% demo script to explain model errors and correlations
% for inflation and randomization (1-dim case)

n = 1000 ; r = 0.6 ;

x = randn(n, 1) ;                          % predictor
e = randn(n, 1) ;                          % noise
y_obs = r*x + sqrt(1-r^2)*e ;              % true predictand

a = x \ y_obs ;                            % regression coefficient (--> r)

y_rgr = x * a ;                            % normal regression
y_ifl =  y_rgr * std(y_obs)/std(y_rgr) ;   % inflation (y_ifl = x in this simple case)
y_rnd =  y_rgr + sqrt(1-a^2)*randn(n,1) ;  % randomization

v_rgr = var(y_rgr) ;                       % equals r^2
v_ifl = var(y_ifl) ;                       % equals 1
v_rnd = var(y_rnd) ;                       % equals 1
					   
r_rgr = corr(y_rgr, y_obs) ;               % equals r
r_ifl = corr(y_ifl, y_obs) ;               % equals r
r_rnd = corr(y_rnd, y_obs) ;               % equals r^2

ev_rgr = 1 - var(y_rgr-y_obs)/var(y_obs) ; % equals r^2
ev_ifl = 1 - var(y_ifl-y_obs)/var(y_obs) ; % equals 1 - 2(1-r)
ev_rnd = 1 - var(y_rnd-y_obs)/var(y_obs) ; % equals 1 - 2(1-r^2)

fprintf('rho = %3.2f\ttrue\test.\ttrue\test.\ttrue\test.\n', r) ;
fprintf('\t\tRGR\t\tIFL\t\tRND\n') ;
fprintf([repmat('-', 1, 60) '\n']) ;
fmt = '\t% 3.2f\t% 3.2f\t% 3.2f\t% 3.2f\t% 3.2f\t% 3.2f' ;
fprintf(['variance:' fmt '\n'], [r^2 v_rgr 1 v_ifl 1 v_rnd]) ;
fprintf(['corr. obs:' fmt '\n'], [r r_rgr r r_ifl r^2 r_rnd]) ;
fprintf(['expl. var.:' fmt '\n'], [r^2 ev_rgr 1 - 2*(1-r) ev_ifl 1 - 2*(1-r^2) ev_rnd]) ;

%% display table
if ~exist('uitable', 'builtin'), exit ; end
cnames = {'true' 'est.' 'true' 'est.' 'true' 'est.'} ;
rnames = {'' 'variance','corr. obs','expl. var.'} ;
dat(1,:) = {'RGR' '' 'IFL' '' 'RND' ''} ;
dat(2,:) = mat2cell([r^2 v_rgr 1 v_ifl 1 v_rnd], 1, ones(6,1)) ;
dat(3,:) = mat2cell([r r_rgr r r_ifl r^2 r_rnd], 1, ones(6,1)) ;
dat(4,:) = mat2cell([r^2 ev_rgr 1 - 2*(1-r) ev_ifl 1 - 2*(1-r^2) ev_rnd], 1, ones(6,1)) ;

f = figure('Position',[200 200 650 130], 'Name', sprintf('RGR-IFL-RND comparison for r = %3.2f', r));
t = uitable('Parent',f,'Data',dat,'ColumnName',cnames,'RowName',rnames);
set(t,'Position',[20 20 600 100], 'ColumnFormat', {'bank' 'bank' 'bank' 'bank' 'bank' 'bank'});
